import numpy as np
import vedo
import vtk
import openseespy.opensees as ops

# 如下为OpenSees中的单元类识别号，用来判断单元类型
ELE_TAG_Truss = (12, 13, 14, 15, 16, 17, 18, 169)  # 169 CatenaryCable
ELE_TAG_Link = (19, 20, 21, 22, 23, 24, 25, 26, 260, 27,  # zeroLength
                86,  # 86为 twoNodeLink
                84, 85, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104,
                105, 106, 107, 108, 109, 130, 131, 132, 147, 148, 149, 150, 151, 152, 153,
                158, 159, 160, 161, 165, 166)  # Bearing
ELE_TAG_Plane = (31, 32, 33, 40, 47, 50, 52, 53, 54, 55, 59, 61, 119, 120, 126, 134,
                 156, 157, 167, 168, 173, 174, 175, 203, 204, 207, 208, 209)
ELE_TAG_Joint = (71, 72, 81, 8181, 82, 83,)
ELE_TAG_Tetrahedron = (179,)
ELE_TAG_Brick = (36, 37, 38, 39, 41, 42, 43, 44, 45, 46, 48, 49, 51, 56, 57, 58,
                 121, 122, 127)


class model_base:
    """
    模型数据获取基类，不要实例化
    """

    def __init__(self):
        pass

    def _get_model_data(self):
        """
        用来获取有限元模型的数据
        :return: None
        """
        # 创造一个虚假的分析对象以方便后面求节点自由度
        ops.wipeAnalysis()
        ops.system('FullGeneral')
        ops.numberer("RCM")
        ops.constraints("Transformation")
        ops.algorithm("Linear")
        ops.integrator('GimmeMCK', 1.0, 1.0, 1.0)
        ops.analysis('Transient')
        _ = ops.analyze(1, 0.0)

        # Get all the node tags
        NodeTags = ops.getNodeTags()
        NodeTags.sort()
        num_node = len(NodeTags)
        # Get all the ele tags
        EleTags = ops.getEleTags()
        EleTags.sort()
        num_ele = len(EleTags)
        # Dict of node coordinates
        Node_Coords = np.zeros((num_node, 3))
        Node_Index = dict()  # key: nodeTag, value: index in Node_Coords
        NodeFix_Coords = []  # nodeTags fix
        model_dims = []

        for i, Tag in enumerate(NodeTags):
            coord = ops.nodeCoord(Tag)
            model_dim = len(coord)
            if model_dim == 1:
                coord.extend([0, 0])
            elif model_dim == 2:
                coord.extend([0])
            model_dims.append(model_dim)
            Node_Coords[i] = np.array(coord)
            Node_Index[Tag] = i
            # dofs
            dofs = ops.nodeDOFs(Tag)
            if dofs.count(-1) == len(dofs):  # all -1
                NodeFix_Coords.append(coord)
        points = Node_Coords
        NodeFix_Coords = np.array(NodeFix_Coords)

        # lines and faces
        # 单元的连接方式，节点数目n-节点号1-节点号2-...-节点号n
        truss_cells = []
        link_cells = []
        lines_cells = []
        all_lines_cells = []
        planes_cells = []
        tetrahedron_cells = []
        brick_cells = []
        ele_midpoints = []  # 单元中点坐标
        for i, ele in enumerate(EleTags):
            ele_nodes = ops.eleNodes(ele)
            # 以下根据单元节点数目来判断单元类型
            if len(ele_nodes) == 2:
                nodeI, nodeJ = ele_nodes
                idxI, idxJ = Node_Index[nodeI], Node_Index[nodeJ]
                all_lines_cells.append([idxI, idxJ])
                if ops.getEleClassTags(ele)[0] in ELE_TAG_Truss:
                    truss_cells.append([idxI, idxJ])
                elif ops.getEleClassTags(ele)[0] in ELE_TAG_Link:
                    link_cells.append([idxI, idxJ])
                else:
                    lines_cells.append([idxI, idxJ])
                ele_midpoints.append((Node_Coords[idxI] + Node_Coords[idxJ]) / 2)
            elif len(ele_nodes) == 3:
                nodeI, nodeJ, nodeK = ops.eleNodes(ele)
                idxI, idxJ, idxK = Node_Index[nodeI], Node_Index[nodeJ], Node_Index[nodeK]
                planes_cells.append([idxI, idxJ, idxK])
                ele_midpoints.append((Node_Coords[idxI] + Node_Coords[idxJ] + Node_Coords[idxK]) / 3)
            elif len(ele_nodes) == 4:
                nodeI, nodeJ, nodeK, nodeL = ops.eleNodes(ele)
                idxI, idxJ, idxK, idxL = Node_Index[nodeI], Node_Index[nodeJ], Node_Index[nodeK], Node_Index[nodeL]
                if ops.getEleClassTags(ele)[0] in ELE_TAG_Tetrahedron:  # 四面体
                    tetrahedron_cells.append([idxI, idxJ, idxK])
                    tetrahedron_cells.append([idxI, idxJ, idxL])
                    tetrahedron_cells.append([idxI, idxK, idxL])
                    tetrahedron_cells.append([idxJ, idxK, idxL])
                else:
                    planes_cells.append([idxI, idxJ, idxK, idxL])
                ele_midpoints.append(
                    (Node_Coords[idxI] + Node_Coords[idxJ] + Node_Coords[idxK] + Node_Coords[idxL]) / 4)
            elif len(ele_nodes) == 8:
                node1, node2, node3, node4, node5, node6, node7, node8 = ops.eleNodes(ele)
                idx1, idx2, idx3, idx4, idx5, idx6, idx7, idx8 = (
                    Node_Index[node1], Node_Index[node2], Node_Index[node3], Node_Index[node4],
                    Node_Index[node5], Node_Index[node6], Node_Index[node7], Node_Index[node8])
                brick_cells.append([idx1, idx2, idx3, idx4, idx5, idx6, idx7, idx8])
                ele_midpoints.append((Node_Coords[idx1] + Node_Coords[idx2] + Node_Coords[idx3] +
                                      Node_Coords[idx4] + Node_Coords[idx5] + Node_Coords[idx6] +
                                      Node_Coords[idx7] + Node_Coords[idx8]) / 8)
        ele_midpoints = np.array(ele_midpoints)

        # 自动根据模型节点坐标判断坐标轴边框的位置
        minNode = np.min(points, axis=0)
        maxNode = np.max(points, axis=0)
        space = (maxNode - minNode) / 15
        minNode = minNode - 2 * space
        maxNode = maxNode + 2 * space
        bounds = [minNode[0], maxNode[0], minNode[1], maxNode[1], minNode[2], maxNode[2]]
        max_bound = np.max(maxNode - minNode)

        # 有限元数据，包括点，线，面，实体的点连接
        self.model_data_finished = True
        self.points = points
        self.NodeFix_Coords = NodeFix_Coords
        self.truss_cells = truss_cells
        self.link_cells = link_cells
        self.lines_cells = lines_cells
        self.all_lines_cells = all_lines_cells
        self.planes_cells = planes_cells
        self.tetrahedron_cells = tetrahedron_cells
        self.brick_cells = brick_cells
        self.ele_midpoints = ele_midpoints
        self.bounds = bounds
        self.max_bound = max_bound
        self.num_node = num_node
        self.num_ele = num_ele
        self.NodeTags = NodeTags
        self.EleTags = EleTags
        self.model_dims = model_dims

    def _eigen_analyze(self, modeTag, alpha):
        """
        本函数用来计算特征值，并保存相应数据
        :param modeTag: 标量，模态号
        :param alpha: 标量，缩放因子
        :return: None
        """
        mode_tag = modeTag
        ops.wipeAnalysis()
        if mode_tag == 1:
            eigenValues = ops.eigen(2)[:1]
        else:
            eigenValues = ops.eigen(mode_tag)
        omega = np.sqrt(eigenValues)
        f = omega / (2 * np.pi)

        if not self.model_data_finished:
            self._get_model_data()

        node_deform = {}
        eigen_vector = {}
        x_eigen = np.zeros(self.num_node)
        y_eigen = np.zeros(self.num_node)
        z_eigen = np.zeros(self.num_node)
        deform_points = np.zeros((self.num_node, 3))

        if alpha is None:
            for i, Tag in enumerate(self.NodeTags):
                coord = ops.nodeCoord(Tag)
                eigen = ops.nodeEigenvector(Tag, mode_tag)
                if len(coord) == 1:
                    coord.extend([0, 0])
                    eigen.extend([0, 0])
                elif len(coord) == 2:
                    coord.extend([0])
                    eigen = eigen[:2]
                    eigen.extend([0])
                else:
                    eigen = eigen[:3]
                x_eigen[i], y_eigen[i], z_eigen[i] = eigen
            eigen_max = np.max(np.sqrt(x_eigen ** 2 + y_eigen ** 2 + z_eigen ** 2))
            alpha = self.max_bound / 10 / eigen_max

        for i, Tag in enumerate(self.NodeTags):
            coord = ops.nodeCoord(Tag)
            eigen = ops.nodeEigenvector(Tag, mode_tag)
            if len(coord) == 1:
                coord.extend([0, 0])
                eigen.extend([0, 0])
            elif len(coord) == 2:
                coord.extend([0])
                eigen = eigen[:2]
                eigen.extend([0])
            else:
                eigen = eigen[:3]
            x_eigen[i], y_eigen[i], z_eigen[i] = eigen
            node_deform[Tag] = np.array([eigen[ii] * alpha + coord[ii] for ii in range(3)])
            eigen_vector[Tag] = np.array(eigen)
            deform_points[i] = np.array([eigen[ii] * alpha + coord[ii] for ii in range(3)])

        eleEigen_midpoints = []  # 变形后的单元中点坐标
        for ii, ele in enumerate(self.EleTags):
            EleNodeList = ops.eleNodes(ele)
            # 以下根据单元节点数目来判断单元类型
            if len(EleNodeList) == 2:
                nodeI, nodeJ = ops.eleNodes(ele)
                eleEigen_midpoints.append((node_deform[nodeI] + node_deform[nodeJ]) / 2)
            elif len(EleNodeList) == 3:
                nodeI, nodeJ, nodeK = ops.eleNodes(ele)
                eleEigen_midpoints.append((node_deform[nodeI] + node_deform[nodeJ] + node_deform[nodeK]) / 3)
            elif len(EleNodeList) == 4:
                nodeI, nodeJ, nodeK, nodeL = ops.eleNodes(ele)
                eleEigen_midpoints.append(
                    (node_deform[nodeI] + node_deform[nodeJ] + node_deform[nodeK] + node_deform[nodeL]) / 4)
            elif len(EleNodeList) == 8:
                node1, node2, node3, node4, node5, node6, node7, node8 = ops.eleNodes(ele)
                eleEigen_midpoints.append(
                    (node_deform[node1] + node_deform[node2] + node_deform[node3] + node_deform[node4]
                     + node_deform[node5] + node_deform[node6] + node_deform[node7] + node_deform[node8]) / 8)
        eleEigen_midpoints = np.array(eleEigen_midpoints)

        self.mode_tag = mode_tag
        self.frequency = f
        self.eigen_points = deform_points
        self.eleEigen_midpoints = eleEigen_midpoints
        self.xyz_eigen = np.sqrt(x_eigen ** 2 + y_eigen ** 2 + z_eigen ** 2)


class OpenSeesVis(model_base):
    def __init__(self, point_size=3, line_width=None, renderLinesAsTubes=True,
                 colors_dict=None, notebook=False):
        """
        point_size: 点尺寸
        line_width: 如果 renderLinesAsTubes为 True 则为 tube半径，默认根据模型边界距离的 1/200
                    如果 renderLinesAsTubes为 False 则为 线宽，默认为4
        renderLinesAsTubes: 将线单元渲染成tube, 此时线宽line_width为tube半径，应根据模型尺寸设置,
                            如果线单元过多，渲染速度可能较慢
        colors_dict: 颜色字典, 默认 dict(point='black', line='#0165fc', face='#06c2ac', solid='#f48924',
                      truss="#7552cc", link="#00c16e")
        notebook: bool, 是否嵌入notebook中显示, 目前有bug，不要打开
        """
        super().__init__()
        self.point_size = point_size  # 节点尺寸
        self.line_width = line_width  # 线宽
        self.renderLinesAsTubes = renderLinesAsTubes
        if notebook:
            self.backend = 'k3d'
        else:
            self.backend = None

        self.model_data_finished = False  # 是否已提取模型数据
        self.mode_tag = None

        # 初始化颜色字典
        colors = dict(point='black', line='#0165fc', face='#06c2ac', solid='#f48924',
                      truss="#7552cc", link="#00c16e")
        if colors_dict is not None:
            colors.update(colors_dict)
        self.color_point = colors['point']
        self.color_line = colors['line']
        self.color_face = colors['face']
        self.color_solid = colors['solid']
        self.color_truss = colors['truss']
        self.color_link = colors['link']

    def model_vis(self, node_label=True, ele_label=True, outline=True,
                  show_fix_point=False, fix_point_size=None):
        """
        本函数用来可视化模型，
        :param node_label: bool，是否显示节点编号
        :param ele_label:  bool，是否显示单元编号
        :param show_fix_point:  bool，是否显示约束点编号， 所有自由度被约束的点被显示
        :param fix_point_size:  float, 约束点尺寸，默认取 point_size 的1.5倍。
        :param outline: bool, 是否显示坐标轴
        :return: None
        """
        if not self.model_data_finished:
            self._get_model_data()
        if fix_point_size is None:
            fix_point_size = self.point_size * 1.5

        if self.line_width is None:
            if self.renderLinesAsTubes:
                self.line_width = self.max_bound / 200
            else:
                self.line_width = 3

        plotter = vedo.Plotter(backend=self.backend, title='OpenSees Viewer')
        # plotter.background_color = "white"

        # all points
        point_plot = vedo.Points(self.points)
        point_plot.color(c=self.color_point, alpha=None)
        point_plot.pointSize(value=self.point_size)
        point_plot.renderPointsAsSpheres(value=True)
        plotter.add(point_plot)
        # show fix
        if show_fix_point:
            fix_point_plot = vedo.Points(self.NodeFix_Coords)
            fix_point_plot.color(c=self.color_point, alpha=None)
            fix_point_plot.pointSize(value=fix_point_size)
            fix_point_plot.renderPointsAsSpheres(value=False)
            plotter.add(fix_point_plot)

        if len(self.truss_cells) > 0:
            if self.renderLinesAsTubes:
                for cell in self.truss_cells:
                    point1 = self.points[cell[0]]
                    point2 = self.points[cell[1]]
                    tube = vedo.Tube([point1, point2], r=self.line_width, c=self.color_truss)
                    plotter.add(tube)
            else:
                startPoints = []
                endPoints = []
                for cell in self.truss_cells:
                    startPoints.append(self.points[cell[0]])
                    endPoints.append(self.points[cell[1]])
                truss_plot = vedo.Lines(startPoints, endPoints=endPoints,
                                        c=self.color_truss, alpha=1, lw=self.line_width)
                truss_plot.smooth()
                plotter.add(truss_plot)

        if len(self.link_cells) > 0:
            startPoints = []
            endPoints = []
            for cell in self.link_cells:
                startPoints.append(self.points[cell[0]])
                endPoints.append(self.points[cell[1]])
            link_plot = vedo.Lines(startPoints, endPoints=endPoints,
                                   c=self.color_link, alpha=1, lw=1)
            plotter.add(link_plot)

        if len(self.lines_cells) > 0:
            if self.renderLinesAsTubes:
                for cell in self.lines_cells:
                    point1 = self.points[cell[0]]
                    point2 = self.points[cell[1]]
                    tube = vedo.Tube([point1, point2], r=self.line_width, c=self.color_line)
                    plotter.add(tube)
            else:
                startPoints = []
                endPoints = []
                for cell in self.lines_cells:
                    startPoints.append(self.points[cell[0]])
                    endPoints.append(self.points[cell[1]])
                line_plot = vedo.Lines(startPoints, endPoints=endPoints,
                                       c=self.color_line, alpha=1, lw=self.line_width)
                plotter.add(line_plot)

        if len(self.planes_cells) > 0:
            verts = self.points
            faces = self.planes_cells
            face_plot = vedo.Mesh([verts, faces], c=self.color_face, alpha=0.6)
            face_plot.lineColor(lc='black')
            face_plot.lineWidth(lw=self.line_width / 5)
            # face_plot.wireframe(value=True)
            plotter.add(face_plot)

        if len(self.tetrahedron_cells) > 0:
            verts = self.points
            cells = self.tetrahedron_cells
            tet_plot = vedo.Mesh([verts, cells], c=self.color_solid, alpha=0.6)
            tet_plot.lineColor(lc='black')
            tet_plot.lineWidth(lw=self.line_width / 5)
            plotter.add(tet_plot)

        if len(self.brick_cells) > 0:
            celltypes = np.empty(len(self.brick_cells), dtype=np.uint8)
            celltypes[:] = vtk.VTK_HEXAHEDRON
            # solid_cells = np.array(self.brick_cells).ravel()
            solid_Data = vedo.UGrid([self.points, self.brick_cells, celltypes])
            solid_Data.alpha(0.6)
            solid_Data.color(self.color_solid)
            solid_Data.lineColor(lc='black')
            solid_Data.lineWidth(lw=self.line_width / 5)
            plotter.add(solid_Data)

        if node_label:
            node_labels = [str(i) for i in self.NodeTags]
            labs = point_plot.labels(node_labels, font='VTK').c('black')
            plotter.add(labs)
        if ele_label:
            ele_labels = [str(i) for i in self.EleTags]
            midpoint_plot = vedo.Points(self.ele_midpoints)
            midpoint_plot.pointSize(value=0)
            labs = midpoint_plot.labels(ele_labels, font='VTK').c('red')
            plotter.add(labs)

        txt = vedo.Text2D(txt='Num. of Node: {0} \n Num. of Ele:{1}'.format(self.num_node, self.num_ele),
                          pos='top-right')
        plotter.add(txt)
        if outline:
            # box frame
            b = vedo.Box(size=self.bounds).alpha(0.1)
            axs = vedo.Axes(b, c='k')  # returns Assembly object
            plotter.add(axs)
        if np.max(self.model_dims) <= 2:
            plotter.show(axes=4)
        else:
            plotter.show(axes=4, viewup='z', zoom=1)

    def eigen_vis(self, modeTag=1, alpha=None, outline=False,
                  show_origin=True, color_map='jet'):
        """
        :param modeTag: int, 模态号
        :param alpha: float, 缩放因子，默认值根据最大变形取模型边界的1/5
        :param outline: bool, 是否显示坐标轴
        :param show_origin: bool, 是否显示原始形状
        :param color_map: 云图类型
        """
        self._eigen_analyze(modeTag, alpha)

        plotter = vedo.Plotter(backend=self.backend, title='OpenSees Viewer')

        # all points
        point_plot = vedo.Points(self.eigen_points)
        point_plot.color(c=self.color_point, alpha=None)
        point_plot.pointSize(value=self.point_size)
        point_plot.renderPointsAsSpheres(value=True)
        point_plot.celldata["mycellscalars"] = self.xyz_eigen  # add an array of scalars to mesh
        point_plot.cmap(color_map, self.xyz_eigen, on='points')
        point_plot.addScalarBar("Eigenvalues")
        plotter.add(point_plot)

        if len(self.all_lines_cells) > 0:
            if show_origin:
                startPoints = []
                endPoints = []
                for cell in self.all_lines_cells:
                    startPoints.append(self.points[cell[0]])
                    endPoints.append(self.points[cell[1]])
                line_plot_origin = vedo.Lines(startPoints, endPoints=endPoints, alpha=1,
                                              lw=self.line_width / 3)
                line_plot_origin.c('gray')
                plotter.add(line_plot_origin)
            # cmap
            if self.renderLinesAsTubes:
                for cell in self.all_lines_cells:
                    point1 = self.eigen_points[cell[0]]
                    point2 = self.eigen_points[cell[1]]
                    scals1 = [self.xyz_eigen[cell[0]]] * 12
                    scals2 = [self.xyz_eigen[cell[1]]] * 12
                    scals = scals1 + scals2 + scals1 + scals2
                    tube = vedo.Tube([point1, point2], r=self.line_width)
                    vmin = np.min(self.xyz_eigen)
                    vmax = np.max(self.xyz_eigen)
                    tube.cmap(color_map, scals, on='points', vmin=vmin, vmax=vmax)
                    plotter.add(tube)
            else:
                # cmap
                startPoints = []
                endPoints = []
                scals = []
                for cell in self.all_lines_cells:
                    startPoints.append(self.eigen_points[cell[0]])
                    endPoints.append(self.eigen_points[cell[1]])
                    scals.append(self.xyz_eigen[cell[0]])
                    scals.append(self.xyz_eigen[cell[1]])
                lines_plot = vedo.Lines(startPoints, endPoints=endPoints, alpha=1, lw=self.line_width)
                lines_plot.cmap(color_map, scals, on='points')
                plotter.add(lines_plot)

        if len(self.planes_cells) > 0:
            if show_origin:
                face_plot_origin = vedo.Mesh([self.points, self.planes_cells], c='gray', alpha=1)
                face_plot_origin.lineWidth(lw=self.line_width / 10)
                face_plot_origin.wireframe(value=True)
                plotter.add(face_plot_origin)
            verts = self.eigen_points
            faces = self.planes_cells
            face_plot = vedo.Mesh([verts, faces], c=self.color_face, alpha=1)
            face_plot.lineColor(lc='black')
            face_plot.lineWidth(lw=self.line_width / 10)
            face_plot.cmap(color_map, self.xyz_eigen, on='points')
            # face_plot.wireframe(value=True)
            plotter.add(face_plot)

        if len(self.tetrahedron_cells) > 0:
            if show_origin:
                verts = self.points
                cells = self.tetrahedron_cells
                tet_Origin = vedo.Mesh([verts, cells], alpha=1)
                tet_Origin.lineWidth(lw=self.line_width / 10)
                tet_Origin.color('gray')
                tet_Origin.wireframe(value=True)
                plotter.add(tet_Origin)
            verts = self.eigen_points
            cells = self.tetrahedron_cells
            tet_plot = vedo.Mesh([verts, cells], c=self.color_solid, alpha=1)
            tet_plot.cmap(color_map, self.xyz_eigen, on='points')
            tet_plot.lineColor(lc='black')
            tet_plot.lineWidth(lw=self.line_width / 5)
            plotter.add(tet_plot)

        if len(self.brick_cells) > 0:
            celltypes = np.empty(len(self.brick_cells), dtype=np.uint8)
            celltypes[:] = vtk.VTK_HEXAHEDRON
            # solid_cells = np.array(self.brick_cells).ravel()
            if show_origin:
                solid_Origin = vedo.UGrid([self.points, self.brick_cells, celltypes])
                solid_Origin.lineWidth(lw=self.line_width / 10)
                solid_Origin.color('gray')
                solid_Origin.wireframe(value=True)
                plotter.add(solid_Origin)
            solid_Data = vedo.UGrid([self.eigen_points, self.brick_cells, celltypes])
            solid_Data.cmap(color_map, self.xyz_eigen, on='points')
            plotter.add(solid_Data)
        txt = vedo.Text2D(txt='Mode {}\nf = {:.3f} Hz\nT = {:.3f} s'.format(self.mode_tag,
                                                                            self.frequency[-1],
                                                                            1 / self.frequency[-1]),
                          pos='top-right')
        plotter.add(txt)
        if outline:
            # box frame
            b = vedo.Box(size=self.bounds).alpha(0.1)
            axs = vedo.Axes(b, c='k')  # returns Assembly object
            plotter.add(axs)
        if np.max(self.model_dims) <= 2:
            plotter.show(axes=4)
        else:
            plotter.show(axes=4, viewup='z', zoom=1)


if __name__ == "__main__":

    def model():

        ops.wipe()

        ops.model('basic', '-ndm', 3, '-ndf', 6)

        ###################################
        # Define Material
        ###################################

        # Define PSUMAT and convert it to plane stress material
        ops.nDMaterial('PlaneStressUserMaterial', 1, 40, 7, 20.7e6, 2.07e6, -4.14e6, -0.002, -0.01, 0.001, 0.3)
        ops.nDMaterial('PlateFromPlaneStress', 4, 1, 1.25e10)

        # Define material for rebar
        ops.uniaxialMaterial('Steel02', 7, 379e6, 202.7e9, 0.01, 18.5, 0.925, 0.15)
        ops.uniaxialMaterial('Steel02', 8, 392e6, 200.6e9, 0.01, 18.5, 0.925, 0.15)

        # Convert rebar material to plane stress/plate rebar
        # Angle 0 is for vertical rebar and 90 is for horizontal rebar
        ops.nDMaterial('PlateRebar', 9, 7, 90.0)
        ops.nDMaterial('PlateRebar', 10, 8, 90.0)
        ops.nDMaterial('PlateRebar', 11, 8, 0.0)

        # Define LayeredShell sections. Section 1 is used for the special boundary elements and section 2 is used for
        # the unconfined interior wall portion
        ops.section('LayeredShell', 1, 10, 4, 0.0125, 11, 0.0002403, 11, 0.0003676, 4, 0.024696, 4, 0.024696, 4,
                    0.024696, 4,
                    0.024696, 11, 0.0003676, 11, 0.0002403, 4, 0.0125)
        ops.section('LayeredShell', 2, 8, 4, 0.0125, 11, 0.0002403, 10, 0.0002356, 4, 0.0495241, 4, 0.0495241, 10,
                    0.0002356,
                    11, 0.0002403, 4, 0.0125)

        # ##################
        # NODES
        # ##################
        # define nodes
        for i in range(11):
            ops.node(16 * i + 1, 0, 0, 0.2 * i)
            ops.node(16 * i + 2, 0.2, 0, 0.2 * i)
            ops.node(16 * i + 3, 0.5, 0, 0.2 * i)
            ops.node(16 * i + 4, 0.8, 0, 0.2 * i)
            ops.node(16 * i + 5, 1.0, 0, 0.2 * i)

            ops.node(16 * i + 6, 0, 1.0, 0.2 * i)
            ops.node(16 * i + 7, 0.2, 1, 0.2 * i)
            ops.node(16 * i + 8, 0.5, 1, 0.2 * i)
            ops.node(16 * i + 9, 0.8, 1, 0.2 * i)
            ops.node(16 * i + 10, 1.0, 1, 0.2 * i)

            ops.node(16 * i + 11, 0, 0.2, 0.2 * i)
            ops.node(16 * i + 12, 0, 0.5, 0.2 * i)
            ops.node(16 * i + 13, 0, 0.8, 0.2 * i)
            ops.node(16 * i + 14, 1, 0.2, 0.2 * i)
            ops.node(16 * i + 15, 1, 0.5, 0.2 * i)
            ops.node(16 * i + 16, 1, 0.8, 0.2 * i)

        ##########################
        # ELEMENTS
        ##########################

        ShellType = "ShellNLDKGQ"
        # ShellType = "ShellMITC4"
        for i in range(10):
            ops.element(ShellType, i + 1, 16 * i + 1, 16 * i + 2, 16 * (i + 1) + 2, 16 * (i + 1) + 1, 1)
            ops.element(ShellType, i + 11, 16 * i + 2, 16 * i + 3, 16 * (i + 1) + 3, 16 * (i + 1) + 2, 2)
            ops.element(ShellType, i + 21, 16 * i + 3, 16 * i + 4, 16 * (i + 1) + 4, 16 * (i + 1) + 3, 2)
            ops.element(ShellType, i + 31, 16 * i + 4, 16 * i + 5, 16 * (i + 1) + 5, 16 * (i + 1) + 4, 1)

            ops.element(ShellType, i + 41, 16 * i + 6, 16 * i + 7, 16 * (i + 1) + 7, 16 * (i + 1) + 6, 1)
            ops.element(ShellType, i + 51, 16 * i + 7, 16 * i + 8, 16 * (i + 1) + 8, 16 * (i + 1) + 7, 2)
            ops.element(ShellType, i + 61, 16 * i + 8, 16 * i + 9, 16 * (i + 1) + 9, 16 * (i + 1) + 8, 2)
            ops.element(ShellType, i + 71, 16 * i + 9, 16 * i + 10, 16 * (i + 1) + 10, 16 * (i + 1) + 9, 1)

            ops.element(ShellType, i + 81, 16 * i + 1, 16 * i + 17, 16 * i + 27, 16 * i + 11, 1)
            ops.element(ShellType, i + 91, 16 * i + 11, 16 * i + 27, 16 * i + 28, 16 * i + 12, 2)
            ops.element(ShellType, i + 101, 16 * i + 12, 16 * i + 28, 16 * i + 29, 16 * i + 13, 2)
            ops.element(ShellType, i + 111, 16 * i + 13, 16 * i + 29, 16 * i + 22, 16 * i + 6, 1)

            ops.element(ShellType, i + 121, 16 * i + 5, 16 * i + 14, 16 * i + 30, 16 * i + 21, 1)
            ops.element(ShellType, i + 131, 16 * i + 14, 16 * i + 15, 16 * i + 31, 16 * i + 30, 2)
            ops.element(ShellType, i + 141, 16 * i + 15, 16 * i + 16, 16 * i + 32, 16 * i + 31, 2)
            ops.element(ShellType, i + 151, 16 * i + 16, 16 * i + 10, 16 * i + 26, 16 * i + 32, 1)

        # P-delta columns

        for i in range(10):
            ops.element('truss', 500 + i + 1, 16 * i + 1, 16 * (i + 1) + 1, 223.53e-6, 7)
            ops.element('truss', 500 + i + 11, 16 * i + 2, 16 * (i + 1) + 2, 223.53e-6, 7)
            ops.element('truss', 500 + i + 21, 16 * i + 4, 16 * (i + 1) + 4, 223.53e-6, 7)
            ops.element('truss', 500 + i + 31, 16 * i + 5, 16 * (i + 1) + 5, 223.53e-6, 7)

            ops.element('truss', 500 + i + 41, 16 * i + 6, 16 * (i + 1) + 6, 223.53e-6, 7)
            ops.element('truss', 500 + i + 51, 16 * i + 7, 16 * (i + 1) + 7, 223.53e-6, 7)
            ops.element('truss', 500 + i + 61, 16 * i + 9, 16 * (i + 1) + 9, 223.53e-6, 7)
            ops.element('truss', 500 + i + 71, 16 * i + 10, 16 * (i + 1) + 10, 223.53e-6, 7)

            ops.element('truss', 500 + i + 81, 16 * i + 11, 16 * (i + 1) + 11, 223.53e-6, 7)
            ops.element('truss', 500 + i + 91, 16 * i + 13, 16 * (i + 1) + 13, 223.53e-6, 7)

            ops.element('truss', 500 + i + 101, 16 * i + 14, 16 * (i + 1) + 14, 223.53e-6, 7)
            ops.element('truss', 500 + i + 111, 16 * i + 16, 16 * (i + 1) + 16, 223.53e-6, 7)

        for i in range(176):
            m = 246000 * 4 / 9.8 / 176
            ops.mass(i + 1, m, m, m, 0, 0, 0)
        # Fix all bottom nodes
        ops.fixZ(0, 1, 1, 1, 1, 1, 1)


    model()

    A = OpenSeesVis(point_size=5, line_width=None)
    A.model_vis()
    A.eigen_vis(modeTag=1, alpha=None, show_origin=True)
    # A.static_vis(num_steps=10, alpha=500, node_label=False, ele_label=False, show_origin=False)
    # A.history_gif(r'C:\Users\yan\Desktop\palte.mp4', num_steps=400, dt=0.1, alpha=100)
